STOCHASTIC EXPONENTIAL INTEGRATORS FOR A FINITE 
ELEMENT DISCRETIZATION OF SPDES 

t3 GABRIEL J. LORD AND ANTOINE TAMBUE 

o 

^«,i Abstract. We consider the numerical approximation of general semilinear 

parabolic stochastic partial differential equations (SPDEs) driven by additive 
^^ space-time noise. In contrast to the standard time stepping methods which 

^ uses basic increments of the noise and the approximation of the exponential 

^^ function by a rational fraction, we introduce a new scheme, designed for finite 

' ' elements, finite volumes or finite differences space discretization, similar to the 

QQ schemes in |4] [7] for spectral methods and [T] for finite element methods. We 

fSJ use the projection operator, the smoothing effect of the positive definite self- 

adjoint operator and linear functionals of the noise in Fourier space to obtain 
I I higher order approximations. We consider noise that is white in time and 

■^r either in H^ or if^ in space and give convergence proofs in the mean square 

y-^ l? norm for a diffusion reaction equation and in mean square H^ norm in 

^^ the presence of an advection term. For the exponential integrator we rely on 

r^ computing the exponential of a non-diagonal matrix. In our numerical results 

we use two different efScient techniques: the real fast Leja points and Krylov 
subspace techniques. We present results for a linear reaction diffusion equation 
in two dimensions as well as a nonlinear example of two-dimensional stochastic 
advection diffusion reaction equation motivated from realistic porous media 
flow. 
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^^ 1. Introduction 

in 

ly-! We consider the strong numerical approximation of Ito stochastic partial difTerential 

f~^ equations 

2 (1.1) dX^[AX^F{X.yX))dt + dW 

> in a Hilbert space H = L'^{n), where fi C M'' and t G [0,r], T > and initial 

data X{0) = Xq is given. The linear operator A is the generator of an analytic 
semigroup S{t) := e^^,t > with eigenfunctions e,, i e N'' and F is a nonlinear 

C^ function. The noise can be represented as a series in the eigenfunctions of the 

covariance operator Q and we assume for convenience that Q and A have the same 
eigenfunctions (without loss the generality as the orthogonal projection can be 
used). In which case [51] we have 

(1.2) I^(x,i)= 5Iv^e,(x)ft(t), 

where q^, i G N'' are the eigenvalues of the covariance operator Q. The Pi are 
independent and identically distributed standard Brownian motions. 
The study of numerical solutions of SPDEs is an active research area and there is 
an extensive literature on numerical methods for SPDEs. Recent work by Jentzen 
and co-workers [21 |3J HI [7] uses the Taylor expansion and linear functionals of the 
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noise for Fourier-Galerkin discretizations of (1.1). In these schemes the diagonal- 
ization of the operator A through the discretization plays a key role. Using a linear 
functional of the noise overcomes the order barrier encountered using a standard 
increment of Wiener process [Ij • In [T] the use of linear functionals of the noise is 
extended to finite-element discretizations (where the operator does not diagonalize) 
with a semi-implicit Euler-Maruyama method. In contrast to [1] here we consider 
two exponential based methods for time-stepping as in [3TJ [5D1 131 [71 [12] . We prove 
a strong convergence result for two versions of the scheme with noise that is white 
in time and in H^ and H^ in space that shows that the exponential integrators 
are more accurate that the semi-implicit Euler-Maruyama method. Furthermore 
we have weaker restrictions on the regularity of the initial data and high accuracy 
for linear problems comparing to the scheme in [1] . The cost of the extra accuracy 
though is that to implement these methods we need to compute the exponential 
of a non-diagonal matrix, which is a notorious problem in numerical analysis |35j . 
However, new developments for both Leja points and Krylov subspace techniques 
[33l (341 \Tl\ (391 \3T\ l32] have led to efficient methods for computing matrix expo- 
nentials. Compared to the Fourier-Galerkin methods of [21 [3l [U [7] we gain the 
flexibility of finite element (or finite volume methods) to deal with complex bound- 
ary conditions and we can apply well developed techniques such as upwinding to 
deal with advection. 



We consider two examples of (1.1 1 where A is the second order operator DA and 



D > is the diffusion coefficient. For the first example 

(1.3) dX ^{DAX - \X)dt + dW 

where A is a constant, we can construct an exact solution. Our second example is 
a stochastic advection diffusion reaction equation in a heterogeneous porous media 

(1.4) dX = (dAX -W -{qX)- — ]dt + dW 

where q is the Darcy velocity field jUj. In the linear example we take Neumann 
boundary conditions and for the example from porous media we take a mixed 
Neumann-Dirichlet boundary condition. 

The paper is organised as follows. In Section [2] we present the two numerical 
schemes based on the exponential integrators and our assumptions on ( |1.1[ ). We 
present and comment on our convergence results. In Section [3] we present the 
proofs of our convergence theorems. We conclude in Section [4] by presenting some 
simulations and discuss implementation of these methods. 

2. Numerical scheme and main results 

We start by introducing our notation. We denote by || • || the norm associated to 
the standard inner product (•, •) of the Hilbert space H — L'^{il) and || • ||_ff'"(o) the 
norm of the Sobolev space 7J™(0), for m G K. For a Banach space V we denote by 
L{V) the set of bounded linear mapping from V to V and L'-^\V) the set of bounded 
bilinear mapping from V x V to C. We introduce further spaces and notation below 
as required. 



Consider the stochastic partial differential equation (1.1), under some technical 



assumptions it is well known (see [231 [231 [2] and references therein) that the unique 



STOCHASTIC EXPONENTIAL INTEGRATORS 3 

mild solution is given by 

(2.1) X{t) = Sit)Xo+ [ S{t-s)F{X{s))ds + Oit), 

Jo 

with the stochastic process O given by the stochastic convolution 

(2.2) 0{t)= [ Sit~s)dW{s). 

Jo 

We consider discretization of the spatial domain by a finite element triangulation. 
Let Tfi be a set of disjoint intervals of fl (for d = 1), a triangulation of fi (for 
d = 2) or a set of tetrahedra (for d — 3). Let V^ denote the space of continuous 
functions that are piecewise linear over Th- To discretize in space we introduce 
two projections. Our first projection operator P/,, is the L^(J7) projection onto Vh 
defined for u e i^(r2) by 

(2.3) iPhU,x) = {u,x) yx^Vh. 

We can then define the operator Ah : Vh — > Vh, the discrete analogue of A, by 

(2.4) {Ah^, x) = {Aip, x) ^,X& Vh. 

We denote by Sh{t) :~ e*'^'* the semigroup generated by the operator Ah- The 
second projection P/v, iV G N is the projection onto a finite number of spectral 
modes e,; defined for u E L^(51) by 

N 

w^ := Pnu = ^(ci, u)ei. 

Furthermore we can project the operator A 

An = PnA and SN{t) := e*^". 

We discretize in space using finite elements and project the noise first onto a finite 
number of modes and then onto the finite element space. The semi-discretized 
version of ( |1.1[ ) is to find the process X'^{t) = X^{.,t) e Vh such that 

(2.5) dX'' = {AhX'' + PhF{X''))dt + PhPNdW, X''(0) = P^Xq. 



The mild solution of (2.5) at time i^ — mAt, Ai > is given by 

X^tm) = Sh{t,n)PhXo + / Sh{t,n ~ s)PhF{X'' {s))ds + / Sh{t,n - s)PhdW^{s). 

Jo Jo 

Given the mild solution at the time f „ , we can construct the corresponding solution 
at tjYi-^1 as 

/.At 

X'^iUn+i) = Sh{At)X''{t„,,)+ Sh{At~s)PhF{X''{s + t^))ds 

Jo 

(2.6) + / "^^ Sh{tm+i - s)PhdW^is). 



For our first numerical scheme (SETDl), we use the following approximations 

F{X''{t^ + s)) « F{X\t^)) se[0, Ai], 
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and 

^' Sh{tm+l~s)PhdW^is) « Ph f"^^' SN{trr.+ l-s)dW^{s) 

= PhPn I "^' 5(t„+i - s)dW{s). 

Then we approximate X^ of X{mAt) by 

(2.7) X::,^, = e^'^'^X':, + AtM^tAh)PhF{X::,) 

+ Ph f "^^ e'-'-^+'-'^^^dW^is) 
where 



rAt 



For efficiency to avoid computing two matrix exponentials we can rewrite the scheme 
( pJl as 

xUi = ^^ + Ai^i(AtA;,) (AX:\ + P,F(X^)) +P„ / "^' e(*-+-^)-4«dW-^(5). 

We call this scheme (SETDl). 

Our second numerical method (SETDO) is similar to the one in [2Tl[5ni[I2. It is 
based on approximating the deterministic integral in ( 2.6 ) at the left-hand endpoint 
of each partition and the stochastic integral as follows 

Sh{tm+i - s)PhdW''is) « Ph f '"^' 5jv(i™+i - s)dW''{s) 

= PhPn I "^' S{t^+^ - s)dW{s). 
With this we can define the (SETDO) approximation YJ^I of X{mAt) by 
(2.8) r,:j+i = M^tAh) (rj: + AtPhF{Y:^J) + Pn /*'"^' e(*™+^"^)^«dW^^(s) 
where 

If we project the eigenfunctions of Q onto the eigenfunctions of the linear operator 
A then by a Fourier spectral method the process 



Ok= / e(*'=+i-^)^"dW^^(s) 



is reduced to an Ornstein-Uhlenbeck process in each Fourier mode as in [1] and 
we therefore know the exact variance in each mode. We comment further on the 
implementation in Section |4] We describe now in detail the assumptions that we 
make on the linear operator A, on our finite element discretization, the nonlinear 
term F and the noise dW. 
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Assumption 1. Let the linear operator —A be a self adjoint positive definite opera- 
tor and A generate an analytic semigroup S . Then there exist sequences of positive 
real eigenvalues {A„}„gpjd and an orthonormal basis in H of eigenf unctions {eiljgfjd 
such that the linear operator —A : 2?(— A) G H ^- H is represented as 

-Av = J2 >^de^,v)e^ V w e V{~A) 

where the domain of —A, 'D{—A) = {w G -ff : ^ Af|(ei,t;)| < 00} and inf Ai > . 

,^fid ten" 

Note that for convenience of presentation we take A to be a second order operator 
as this simpHfies notation for the norm equivalence below. Similar result hold, how- 
ever, for higher order operators. We recall some basic properties of the semigroup 
S{t) generated by A that may be found for example in [T51l25| . 

Proposition 2.1 (J^). Let /3 > and < 7 < 1, then there exist C > such 
that 

\\{-AfS{t)\\ < Ct"P for t>0 

\\{-Ay{L-S{t))\\ < Ct^ for t>Q. 

In addition, 

{-AYS{t) = S{t)(-AY on V{{-Af) 

If P > 1 then V{{~Af) C V{{-Ay'). 

We introduce two spaces H and V where H C F that depend on the choice of the 
boundary conditions. For Dirichlet boundary conditions we let 

y = H = Hl{n) = {veH^{n):v^Q on dVi}, 

and for Robin boundary conditions, for which Neumann boundary conditions are 
a particular case, V — H^{il) and 

H = {w e H\n) : dv/dvA + av = Q on dVl) , 

see [TB] for details. Functions in H satisfy the boundary conditions and with M 
in hand we can characterize the domain of the operator (— A)*"/^ and have the 
following norm equivalence [TU [5S] for r = 1 , 2 

We now introduce our assumptions on the spatial domain and finite element space 
Vh- We consider the space of continuous functions that are piecewise linear over 
the triangulation Th with Vh C V . 

Assumption 2. (Regularity of the domain Vt and the space grid) 

We assume that 17 has a smooth boundary or is a convex polygon and that the 
maximal length h of the elements of Th satisfy the usual regularity assumptions on 
the triangulation i.e. for r — 1,2 

(2.9) mi{\\v-x\\ + h\\v{v~x)\\)<Ch^\HH^in), vevnff'in). 

This inequality is sometimes called the Bramble and Hilbert inequality, see [TB] or 
[H]. It follows that 

(2.10) \\PhV-v\\<Ch^\\v\\Hria) yveVnH^n), r = l,2. 
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Assumption 3. (Nonlinearity) 

Let V be a separable Banach space such that 2?((— A)^") C V C H = L'^{il) 
continuously. We assume that there exist a positive constant L > such that the 
Nemytskii operator F satisfies one of the following 

(a) _F : V — > V is twice continuously Frechet differentiable mapping with 

\\F'{v)w\\<L\\w\\, ||F'(«)|U(v)<i, ||F"(«)|L(.,(v)<L 

and 

II(^'H)11l(d((-a)V2)) < L(i + hllx,((-^)i/.)) ^v,wev, ue v({-Af'% 

where (F'{u))* is the adjoint of F'{u) defined by 

{{F'{u))*v,w) = {v,F'{u)w) Vw,w e H = L'^{n). 
As a consequence 

\\F{X)-F{Y)\\<L\\X-Y\\ yX,YeH, 
andyX e H = L^{n) 
\\FiX)\\ < \\Fm + \\FiX) F(0)|| < ||F(0)|| + L\\X\\ < C(||F(0)|| + ||X||). 

(b) F is globally Lipschitz continuous from {H^{fl), ||-||ifi(o)) to [H = L^(il), ||.||) 
then 

\\F{X) - F{Y)\\ <L\\X- r||Hi(o) VX, Y e H\n). 

We assume that the function F is defined in L^(il), although in general F may be 
defined in any Hilbert space. The possible choice of V can be if, H^{n) or the M— 
Banach space of continuous functions from [0,T] to H denoted by C{[0,T],H) if 
d=l. 

We now turn our attention to the noise term and introduce spaces and notation that 
we need to define the Q- Wiener process. Denoting by C{H) the Banach algebra of 
bounded linear operators on H with the usual norm. We recall that an operator 
T e C{H) is Hilbert-Schmidt if 

If we denote the space of Hilbert-Schmidt operator from Q^^^{H) to H by L" := 
HS{Q^/^{H),H)i.e 

I ieN'^ 



the corresponding norm \\.\\i^o by 



Let (p{u>) be a process such that for every sample lu, (p{oj) G L!]- Then we have the 
following equality 

E|| [\dwf= fnwwi.ds^ f nwQ^'^Wnsds, 

Jo Jo Jo 
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using Ito's isometry J4|. We assume sufRcient regularity of the noise for the exis- 
tence of a mild solution and to project the noise into the finite element space V/j. 
To be specific we assume the noise is in either in H^ or H^ in space. 

Assumption 4. (Regularity of the noise) For all t g [0, T] we assume that 
0{t) is an adapted stochastic process to the filtration {J-t)t>o with continuous sample 
paths such that 0{t2) ~ S{t2 — ti)0{ti) , < ti < t2 < T is independent of J-t^ LetV 
he a separable Banach space such that !?((— A)^' ^) C V C H = L^(J7) continuously. 
For some 9 G (0, f /2] 

0(i) e 2?((-A)'-/2) ^ MnH'-in), r = l,2 

E{\\0{t2)~0{t,)\\^v) < C{t2-t^f'. 

Using the equivalence of norms, we have that 

0{t)^V{{-AY'^) \ft^[Q,T]^\\{-AY/^Q^'^\\Hs<oo r = l,2. 

2.1. Main results. Throughout the paper we let N be the number of terms of 
truncated noise and Xjq = {1,2,. ..,7V} and take t^ — mAt G (0,T], where T = 
MAt for m, M G N. We take C to be a constant that may depend on T and 
other parameters but not on At, N or h. We also assume that when initial data 
Xo G I?((-A)T) then EWi-A^XoW^ < cx), / = 2,4 and < 7 < 1. 
Our first result is a strong convergence result in L^ when the non-linearity satisfies 
the Lipschitz condition of Assumption Is] (a) with scheme (SETDl). This is, for 
example, the case of reaction-diffusion SPDEs. 



Theorem 2.2. Suppose that AssumptionsUi pi ^Va) and\M are satisfied with r 
1,2. Let X{tm) be the mild solution of equation (1.1 1 represented by (2.1) and X 



be the numerical approximation through (2.71 (SETDl). Let < 7 < 1 and set 
a = inm(29,j) and let 9 G (0, 1/2] be defined as in AssumptionlA 
IfXoeV{{-AY) then 



{E\\X{t,,)~X:^J'f < C f-'/^ + Ai'^+f inf A, 



-r/2^ 



IfXo G T>{-A) =MnH^{n) then 

{E\\X{t^) - A„\||2)'/' <clh^ + At^' + ( . inf A, 



r/2\ 



Our first result for scheme (SETDO) is a strong convergence result in L^ when the 
non-linearity satisfies the Lipschitz condition of Assumptions^ (a). 

Theorem 2.3. Suppose that AssumptionsYn IM ^Na) andu^ are satis fied with r 



1,2. Let X(tm) be the mild solution of equation (1.1) represented by (2.1) and Y^ 



be the numerical approximation through (2.81 (SETDO). Let < 7 < 1 and set 
a = min(20,7) and let 9 G (0, 1/2] be defined as in Assumption^ 
If Xo G V{{~Ay) then 

{E\\X{tra) - Ytf)^''' < C t^^^'^h"- + At'' + At |ln(Ai)| + i Jnf^ \, 
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ifXo e v{-A) = mnH'^{n) then 

h"^ + At'' + At \HAt)\ + ( mf^ A, 

For convergence in the mean square iJ^(ri) norm where the non-hnearity satisfies 
the Lipschitz condition from L^(i7) norm to H^{D,) ( Assumptions^ (b)) we can 
state results for (SETDl) and (SETDO) together. 

Theorem 2.4. Suppose that Assumptions [71 III IWb), Q] (with r ~ 2) are satisfied 

and F{X) G V with E sup ||i^(X(s))||j:^i(Q) ) < cxd. Let X be the solution mild 
\0<s<T ) 



oj equation {1.1) represented by equation {2.1) and Cj^ be the numerical approxima- 
tions through scheme (pjl) or E^ ( C^ ^X^ for scheme (SETDl) and C,^ = Y^ 




-1/2N 



for scheme (SETDO)). Let < 7 < 1- Then we have the following: 
//XoeP((-A)(i+'')/2) then 

{E\\X{tr,)-Ct\\W)Y'^ < c({t;^'/'h + At-'/')+( inf A 

IfXoeV{{-A)) then 

{E\\X{t„,)-Crm^,,^y/' < c({h + At'/'-^)+( inf A, 

for very small e G (0, 1/2). 

We note that this theorem covers the case of advection-diffusion-reaction SPDEs, 
such as that arising in our example from porous media. 

We remark that if we denote by Nh the numbers of vertices in the finite elements 
mesh then it is well known (see for example [37]) that if A^ > N^ then 

inf A, I <Ch' and ( inf A,| <Ch. 



As a consequence the estimates in Theorem |2.2[ Theorem |2.3| and Theorem |2.4| can 
be expressed in function of h and At only and it is the error from the finite element 
approximation that dominates. If A^ < Nh then it is the error from the projection 
P/v of the noise onto a finite number of modes that dominates. 



From Theorem 2.4 we also get an estimate in the root mean square i^(r2) norm 
in the case that the nonlinear function F satisfies Assumption [s] (b) . We cannot 
do the proof directly in L''{^) due to the Lipschitz condition in Assumption ^ (b) . 



Simulations for Theorem 2.4 will be do in L'{Vt) since the discrete L''{Vi) norm is 
more easy to use in all type of boundary conditions. 

Finally if we compare these theorems to those in ^ for a modified semi-implicit 
Euler-Maruyama method then we see that using the exponential based integra- 
tors we have weaker conditions on the initial data and in particular the scheme 
(SETDl) has better convergence properties. 
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3. Proofs of main results 



3.1. Preparatory results. We start by examining the deterministic linear prob- 
lem. Find u E V such that such that 

(3.1) u' ~ Au given u(0) = v. 

The corresponding semi-discretization in space is : find Uh G Vh such that 

where w^ — PhV. Define the operator 

(3.2) Th{t) := S{t) - Sh{t)Ph = e'^ - e'^'^Ph 
so that u{t) — uii{t) == Th{t)v. 

Lemma 3.1. The following estimates hold on the semi-discrete approximation of 

(3.3) \\u{t)~uh{t)\\ = \\n{t)v\\<Ct-^'^h'\\v\\ tf veH = L\n), 

(3.4) \\u{t)-uhm = ||T,(t)t;|| <C;i2||^||^^^^^^ tfveV{{-A}), 

i3.5)\\u{t)-Uh{t)\\HHn) = \\Thit)v\\Hi(n)<Cht-'/^\\v\\Hnn) ifveV, 
i3.6)\\uit)-Uh{t)\\HHn) = \\Th{t)v\\m(^n)<Ch\\v\\H2(^n) ifveVi^A). 



Proof. The proof of the estimates ( |3.3[ ) and (3.4) can be found in 15J with Dirichlet 
boundary conditions. The same proof can be generalized easily to Robin or mixed 
boundary conditions, incorporating the extra term from the boundary with the 
bilinear form. Estimates (3.3)- (3.6) are the special case of the proof of Theorem 
3.1 in |14j where the nonlinearity is taken to be zero. For our case 

u{t) = S{t)v, 

and we have the following estimates for t e (0, T] 

\Ht)\\H^,^n) < ct-^'-^y^v\\Hi(n) a vem 

\Ht)\\H2^n) < C\\v\\H2in) if veV{-A), 
||"t(t)lk=(o) < C7t-i-(^-i)/2||«||^i(f,) if i;eH 
\\utit)\\m(n) < Ct-'/'-\\v\\H2^n) if veVi-A) 
Using these in the proof of [TH Theorem 3.2] gives the result. D 



s — 


1,2, 


s 


= 0,1, 


s - 


= 0,1. 



We now consider the SPDE ( 1.1 1 



Lemma 3.2. Let X be the mild solution of (1.1) given in {2.1), let Q < ^ < 1 and 

ti,i2e[0,r], ii<i2. 

(i) If Xq e V{{-A)^), ||(-A)"/2gi/2||^^ <■ oo „ji/j Q < a <2 and suppose F 
satisfies Assumption^ (a). Set a = min(7, 1/2, a/2) then 

E\\X{t2) - X{ti)f < Cit2-tif"lE\\Xo\\^^ + E( sup (||F(0)|| + ||X(s)||)^ 



VO<s<T 

Furthermore 

E\\{Xit2)-0{t2))-iXit,)-0{t,))f < C{t2-t,f'^ 0<7<1. 
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(ti) If Xa e 2?((-A)(T+i)/2), ||(-A)1/2q1/2||^^ < ^ and F{X) e H\n) with 



E sup ||F(X(s))||hi(o) <(X)then 
\a<s<T 



^ 2 
2 ^ rff-i J- \f ( ip\\ \^ ||2 



E||X(t2) - X{h)\\' < C{h - hV E\\Xo\\l^+,y^ + E sup \\F{X{smH^m 

Y \0<s<T 

Remark 3.3. Before doing the proof it is important to notice that if Xq G 2?((— A)^) 
with E||(— A)'''Xo)||' < oo, I — 2,4, Assumption 1-4 ensure the existence of the 
unique solution X G 'D{{—Ay') such that 

e( sup \\{~AyX{s))\A <oo. 

\0<s<T J 

In general if Xq G X>((-A)'^) and 0{t) G ^((-A)") i/ien X G P((-A)™"(T'°)) 



E ( sup 

vO<s<T 



ll(-A)™"('''")x(s)|n <oo. 



More information about properties of the solution of the SPDE [1.1) can he found 
in [7]. 

Proof. The first claim of part (i) of the Lemma can be found in [T and so we prove 
the second part of (i). Consider the difference 

{X{t2)+0{t2)-{X{t^)+0{t^))) 

= {S{t2) - 5(ti)) Xo+fJ' S{t2 - s)F{X{s))ds - J ' S{ti - s)F{X{s))ds 
= I + 11 
so that 

mX{t2) + 0{t2)) - (X(ti) + 0(ti))||2 < 2(E||/||2 + E||//||2). 



We estimate each of the terms /, // . For < 7 < 1, using Proposition 2.1 yields 

||/|| = \\s{t^){-A)-^{i-s{t2-t^))(-Ayx4 < c{t2-t,y\\Xo\u. 

Then E||/||2 < C(i2 - ti)2'^E||Xo||^. For the term //, we have 

II = f {S{t2 -s)- S{ti - s))F{X{s))ds + f ' 5(t2 - s)FiX{s))ds 

Jo Jti 

= Ih + Ih. 

We now estimate each term Hi and II2. For E||//2|p boundedness of S gives 



EII//2II' < (j\\\S{t2~s)F{X{s))\\ds 
< C{t2-tifE( sup \\F{X{s))\ 



2 



0<s<T 
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For E||//i|p we have 



Ih = I iSit2-s)-S{h-s))FiXis))ds 

' \s{t2 -s)- S{h - s)) {F{X{s)) - F(X{ti))) ds 



+ f\s{t2-s)-S{h-s))F{X{h))ds 
Jo 

= //11 + //12. 

Using the Lipschitz condition in Assumptions^ (a) with the first claim of (i) yields 



E||//iiir < 



(Sih -s)- Sih - s))E||(X(s) - X{h)\\ds 







< c[{t2-h) ih-sy-'ds 

\2 



< cit2-hr 

Assumption [3] (a) gives 

{E\\Ih2\\Y' < {nFiXih)\\Y'\\ riS{t2-s)-Sih~s))ds\\ 

Jo 



< C\\ / 5(t2 - s) - S{h ~ s)ds\\ 

'0 



Using the fact that S is bounded we find 



EII//12I 



< c\ 

< c\ 
= c\ 

= c\ 



''^' -- ^Sih) I {S{t2-h-s)-S{-s))ds\\ 
Jo 

''\s{t2-h + s)^S{s))ds\\ 
Jo 

S{s)ds- f ' S{s)ds\\ 

't2~ti Jo 

S{s)ds + f ' S{s)ds - / ' Sis)ds\\ 

/ta-ti -Jti Jo 

' ' S{s)ds- [' ' S{s)ds\\ 

Hi Jo 

< C{t2~h). 

Combining the previous estimates ends the proof of the second claim of (i) . 
We now prove part (ii) of the lemma. Consider the difference 

X{t2)-Xih) 

= {S{t2) - S{h)) Xo+(r S{t2 - s)F{X{s))ds - / ' Sih - s)F{X{s))ds 

+ ( f' S{t2 - s)dW{s) - f ' S{ti - s)dW{s) 
= I + 11 + III 
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and then 

Let us estimate the terms /, // and /// and we start with /. If Xq G 'D{{-A)'^''+^^^^) 
using Proposition |2.1| yields 

UWrnm = \\i-A)'/'Sih)il-Sih-h))Xo\\ 

= \\{-Ay/^Sih)il - Sih - ti))(-A)-^/2(-^)^/2Xo|| 
= ||5(ti)(-A)-^/2(I - Sih - ii))(-A)(^+i)/2Xo|| 

Then 

ni\\min<Cih-hr\\Xo\\l+,y 
For the term //, we have 

II = [' {S{t2 - s) - S{h - s))F{X{.s))ds + f ' S{t2 - .s)FiX{.s))ds 

Jo Jti 

= Ih + Ih. 

We now estimate each term above. Using the fact that in 2?((— A)^/^) we have the 
equivalency of norm ||.||hi(si = IK^^)"'" -lli "^6 have 

(3.7) \\S{t)U(mm<\\{-A)^/^S{t)\\ 

where ||5'(i)||L(/fi(s7)) is the norm of the semigroup viewed as a bounded operator 
in iy^(J7). We also have the similar relationship for the operator S{ti) — S{t2) with 
h,t2G [0,T]. 
Using similar inequality as (|3.7|) yields 



E||//i||2,,( = E\\ f\s{t2-s)-S{h~s))F{X{s))ds\\l,^^) 
Jo 

< E ('y ' \\{S{t2 - s) - S{h - s))F{X{s))\\HHn)ds] 

< ( f \\{S{t2 -s)- S{t, - s))\\H^n)ds] E ( sup ||F(X(s))||^i(, 

< (J ' \\{^Ay/\S{t2 - s) - S{h - s))\\ds] E(snpJ\F{X{s))\\H^n)] 



For e € (0, 1) small enough, we have 

< 



< 



[ ' ||(-A)(i-^)/25(ti - s)(-A)(^-i)/2((I - S{t2 - hMlds) E ( sup ||F(X(5))|Ui 

/ Vo<s<T 

C{t2 - hy-' ( f\ti- s)^'-^y'ds) E ( sup ||F(X(s))|Ui(o)') 

\J0 J \0<s<T J 



(O) 



< C(t2-ti)i-'E( sup |lF(X(s))||ai(jj)" 

vO<s<T 
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We also have using Proposition |2.1| 

< E ( f ' \\S{t2 - s)F{X{s))\\HUn)ds' 



< E 



L / 

\\Sit2-s)\\LiHHn))\\FiXis))\\H^n)ds 



ti 



< E 



t2 



{-Ay^\S{h-sm\\FiX{s))\\HHn)d.^- 



< 



[\t2--s)'^'^ds\ Ef sup \\F{X{s))\\Hi(n) 

Jti / \0<s<T 



< C{t2-h)-E[ sup \\FiXi.s))\\H^^n) 
\n<s<T 



Hence, if F{X) e if ^(f^) with E sup \\F{X{s))\\HUn) I < oo, we have 

Vo<s<T / 

E||//||2 < 2(E||//i||2 + E||//2|n < C(i2 - ti)^E ( sup ||F(X(s))|Ui(o) 

\0<s<T 

We also have for the term /// 

/// = / S{t2 ~ s)dW{s) - I S{h - s)dW{s) 
Jq Jo 

{S{t2 -s)- S{h - s)) dW{s) + [ ' Sit2 - s)dWis) 

Jti 



= Ilh+IIh. 
The Ito isometry property yields 



E||///i||?,i 



ffi(n) 



= E|| / ' {S{t2 -s)- Sih - ■5))dW{s)\\j,,^^^ 
Jo 

< r E||(-^)i/2 (^S{t2 -s)- Sih - s)) Q'/Ynsds 
Jo 

E|| (5(<2 ~s)~ S{h - -s)) {-A'/')Q'^'\\lsds. 



Using Proposition 2.1 the fact that S{t) is bounded and \\i—Ay/^Q^/^\\HS < oo 
yields 

E||///i|ll,i(o) < C f'\\{S{t2-.s)-S{h~sm^ds 
Jo 

= C f ' \\i-A)'^'-'^/^S{h- s){-A)-^'-'y^l- Sit2-h))\\^ds 



< C{t2-tiY-' / {h^s)-'+'ds 
Jo 

< c{t2-hY''. 
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with e G (0,1) small enough. Let us estimate E||///2||i/i(o)- The fact that 

|j(_A)i/2Qi/2||^^<^ yields 



< f\\{^Af/^S{t^-s)Q'/YHsds 



< C{t2~h). 



Hence 



niiif < 2(Eii///if + Eii///2f ) < c(t2 - hr. 

Combining the estimates of E||/|p, E||//|p and E||///|p ends the proof. D 

Remark 3.4. 7/7 > 1 and with more regularity of the noise (0{t) G V^^—AY), r > 
l/2j we have 



E\\X{t2)-Xiti)f<C{t2-ti) 



l-e 



for any e S (0, 1). 

We can prove that we can take 9 = 1/2 for V = H^{Vi) or if 0{t) e X>(-A). We 
have 9 7^ 1/2 and close to 1/2 if 0{t) G 'D((—A)^''^). These estimates follow those 
used to estimate II Ii and III2 in the proof of Lemma 3.2. see jT5] . 



3.2. Proof of Theorem 2.2, The proof follows the same basic steps as in [T], 
however here the discrete semigroup is an exponential. As a consequence the es- 
timates are different and the proof here is simpler with fewer terms to estimate. 
Set 

m-l „tk+i 
X{tm) = S{t,n)Xo+Y. S{t,n-s)F{X{s))ds + 0{t^) 

k=0 ■^*'' 

= X{t,n)+0{t„,). 

Recall that by construction 

X^^ = e^'^'^X';^^_,+ e^'^'~'^^'^PhF(X^^_,)ds + Ph e^''--''>^"dW^ {s) 

Jo Jtm-i 

™^l / rtk+i ftk+i \ 

= Sh{tm)PkX„ + II / ^h{t,n - s)PhF{X^)ds + Ph / SN{tm " s)dW'' {s) 

m-l , tk + i \ 

= Sh{t^)PhXa +J2[ '^''(^" - s)PhF{Xli)ds J + PhP^Oit^) 

k=Q ^•^*'= ^ 

^ ^m + PhPNO{tm), 
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where 



m-1 . „t^._|_j 
= Sh{t„^)PkX„ +Y.( ^hiUa - s)PhF{Zt + PhPNO{tk))ds 

fc=0 ^•^*'= 



1 /q 

We now estimate {E\\X{tm) - X'^W^) . We obviously have 

X{tm) — Xj^ — Xt^ + 0{t,n) — XjI^ 

= X{tra) + 0{tm) ~ {Z'^ + PhPNO{t„,)) 

= (X(t„0 - Z,';) + (PAr(O(i„0) - PhPN{0{t„,))) + {Oit„,) - PAr(O(t„0)) 

(3.8) = I + 11 + III. 

Then (E|lX(t„) - X^^f)'^^ < {mrf + (E|1//|P)'/' + (E|1///||2)^/' and we 
estimate each term. Since the first term will require the most work we first estimate 
the other two. 

1 /2 

Let us estimate (E||//|p) . Using the property (2.101 of the projection P^, the 
equivalence |i.|U.-(o) = U-^T^^-W in V{{-AY^^), the Ito isometry and the fact 
that the semigroup is a bounded operator yields 



E||//||2 < Ch^'-mi-AY^^ S{Un-s)dW{s)\\ 

Jq 

< Ch^' f '" \\{-AY^^S{t^ - s)\\lods 

Jo ^ 

< ch'^ r Wi^Ay/'Q'fYHsds. 



2 







Thus, since the noise is in iP' we have (E||//||2)i/2 < ch'' 
For the third term /// 



E||///||2 = E||(I - PN)Oitm)\\' = E||(I - PAr)(-A)--/2(-A)'-/20(i,„)||2, 



and so 



E\\IIlf<\\il-P^)i-A)-^/Ym-Ay/'0{t^)\\'<c{ inf A, 
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We now turn our attention to the first term E||/|p. Using the definition of T^ from 
(3.2) the first term / can be expanded 

™-i j-tk+i 

I = r,,Xo + I] / S{tm ' S)F{X(S)) - Sh{tm - s)PhF{Z^ + PhPNO{tu))ds 

k=0 •^*'=* 

= nXo +J2 ^''(^" ~ s)Pi,{F{X{tk)) - F{Zi: + PhPNO{tu))))ds 

k=0 -'**' 

™-l j-tk+1 
+ Y, / Sh{t-^-s)Ph{F{X{s))-F{X{tu)))ds 

k=0 •'*i-= 
™-l rtk+i 
+ J2 / iS{t^ - s) - Sh{tm - s)Ph)F{X{s))ds 

= h + h + h + h- 
(3.9) 



Then 



(E\\I\ff' < (EWhlff 



For /i, if Xo e P((-yl)T) C i/, equation ([33| of Lemma [sl] gives 

gives 



and if Xq G X>(-yl) =mr\H^{n), equation (3.4) of Lemma 



3.1 



{nhfy^'<ch'(nxo\\H' 



2(0) 



1/2 



If i^ satisfies Assumption^ (a), then using the Lipschitz condition, triangle inequal- 
ity as well as that Sh{t) and P^ are bounded operators, we have 



{nh\n 



2\l/2 



>l/2 



™-l /.tfc + i 

< Cj] / (E||i^(X(ifc))-i^((^fe^ + ^/.i'^0(tfc))f)^''ds 

in — l ptk-^i 

< cY, / {nx{tk)-xtrf''ds. 



./M|2^1/2 

fc=0 •'^'= 

As /3 needs more work let us estimate I^ first. Using the fact P^, S", S*/! are bounded 
with (3.3) of Lemma 3.1 yields 

2\l/2 



{nhry 



m-l /.tfc + i 

fc=0 •^*'= 



'E||T;,(i^-s)F(X(s))|p)'/'ds 



< C/i^ sup (E||F(X(s))| 

0<s<T 



2\l/2 



(t. 



-1/2 



< Ch' 
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Let us estimate (E||/3||^) . We add in and subtract out 0{s) and O(ffc) yields 



fc=0 •'^'= 



= ( E|| Y. I ^h{t,n - s)Ph {F{Xis)) - FiX{tk) + Ois) - OihM ds\ 
\ k=o •^*'= / 

/ m-l „tj^^j \ 

+ I E|l ^ J Shitm - s)Ph (F(X(tfe) + 0{s) - Oih)) - FiXih)))) ds\\^ J 

:= {E\\llff' + E{\\l!\ff". 

Applying the Lipschitz condition in Assumption p[ a), using the fact the semigroup 
is bounded and according to Lemma 3.2 for Xq e 'D{{—A)''), < 7 < 1 we 
therefore have 



1/2 



inmT' < ^E / {mx{s)-o{s))-{xit,)-o{t,))\ 



1 /q 

Let us now estimate E (||/|||^) . The analysis below follows the same steps as in 
[7] although the approximating semigroup Sh is different here. Applying a Taylor 
expansion to F gives 

with 

/|i = E|| ^ / Shit„, - s)P„F'(X(ifc))(0(s) - Sis ~ ifc)0(tfc))ds||2 
V fc=0 •^*'' / 

If = [nJ2 Shitrn-s)PhF'{Xitk))iSis-tk)0{tk)-Oitk))dsf\ 

\ k=0 •^*'= / 

/ m-l .tfc + i .1 \ ^/^ 

If = E||^ / Shit„,-s) Ril - r)drds\\^ ] , 

R := F„i^"(X(tfc)) + r(OGs)-0(tfe))(OGs)-0(tfe),0(s)-0(tfc)). 

Using the fact that 0(^2) - S{t2 - ii)O(ti), < fi < ^2 < 7" is independent of J^^, 
one can show, as in [7], that 

(^3')' = Ee|| r^' Sh{t^ - s)PhF'{X{tk)){0{s) ~ S{s - tk)0{tk))dsf. 
fc=0 •^*'' 
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Therefore as Sh is bounded we have 

r21 

-'3 



< {Y.[ {nSh{tm-s)P,,F\X{tk)){0{s)-S{s-tk)0{tk)Wy'd. 

1/2 
/m-1 / /-tfc+i 

< ciY,( {nPhF'{X{tk)){0{s)-Sis-t,)0{tk)W) ds 

\ L — n V''*fe 



Using Fubini's theorem in integration with Assumption Wl Assumption pTa) and 
Proposition |2.1| yields 



r21 



^"^1 ^tfc+i ^ ^/^ 



< CAti/2 ^ / E\\PhF'{X{tk)){0{s) - S{s - t,)0{tk)Wds 

< CAf/^ Y. / E||(0(s) - S{s - tk)Oitk))fds 



Cm— 1 ^tfc_|_i 2 ^ 

^ / ((E||0(s)-0(tfc)||2)'/%(E||(5(s-tfc)-I)0(tfe))||2)'/') d5 

fe=o •^*'' y 

< cAt'/'iY.l {{^ - hf + is - hY^' {noit,mf') dsj 

Let us estimate /|^. 

t22 

-'3 

'"-I /-ifc+i , .1/2 



1/2 



k=o •^*'= 

$: / ||5,(i™-s)P,(-A)V2|| (E\\{-A)-'/'F'{X{t,)){S{s-tk)-l}0{t,))r) 



k=0 "= 

< i] /" '*' \\S,At^ - s)P,A~AV/^\\ (n(-A)-'^^F'(X(t,,))(S(s - t,) -I)0(UM^Y^^ ds. 

fe=0 

Since P;i(— ^)^/^ = (— A/i)^/^ and 5*^ satisfies the smoothing properties analogous 
to S{t) independently of h (see for example [H]), and in particular 

\\Sh{trn){-AhY^^ - \\{-Ah)^/^Sh{t,n)\\ < ^C^', ^m = mAt > 0, 

we therefore have 



r22 

^3 




< 


k = •'*!' 



< CY, f "^^ (tm - sy^'^ (n{-A)-'/^F'{X{tkmS{s - tk) - I)0(tfe))||2) "^ ds. 

1 n ^ th 
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The usual identification oi H = L'^{VL) to its dual yields 



t22 

-"a 



™-i /-tfc+i 






1/2 



fc=0 ^'^ 



1/2 



E sup |(«, (-A)-i/2F'(X(tfe))((5(s - tk) - l)0(ife))>| 



ds 



\\v\\<l 



where (, ) = (, ) and we change the notation merely to emphasize that H is 
identified to its dual space. The fact that (— A)^^/^ is self-adjoint implies that 
((-A)-i/2^'(X))* = F'(X)*(-A)-i/2 This combined with the fact that H C 
X>((-^)i/2) thus ^((-A)-!/^) c (X>((-A)i/2))* c H* = H continuously and 
Assumption w[a) yields 



™-l rtk+i 



I3 < ^.Z^ / (im ~ i/c+l) 



-1/2 



fe=0 •"^'= 



1/2 



E sup \{F'iX{tk)n~Ay'^\, {S{s - h) - l)0(tfe))| 



ds 






k=0 •'*!' 



-1/2 



1/2 



E sup \\F'{x{tk)n-AryMi\\is{s - tk) -i)o{t,))\u 

\M\<i , 

k=Q •^** 



ds 



.(E(i + iix(t,.)iii)'ii(^(s-t.)-i)o(ife)ii-ir 



1/2 



ds 



fe=0 '^'» 



1/4 



1/4 



ds 



1/4 



ds 



. (e (1 + ||x(i,.)lli)') (e {\\S{s - t,) - l)0{t,)\U) 

m—l 

< cJ2{t'm-tk+iy'^^ 

k=0 

. (1 + {E\\Xit,)\\tf') T"' (e iWSis - t,) - I)O(tfe)ll-i) 

771—1 

\\{-A)-^^'-+'/'HS{s^t,)-l)\\{nO{t,)tf'ds 

< CJ2 (tm - tk+l)'^'^ I '^' ||(-A)l/2--/2(_^)-l(5(, _ tk) - I)||ds 



fe=0 



tk 

771—1 



/c=0 
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Using Proposition 2.1 and the fact that {—Ay/^ ^^^ is bounded as r = 1, 2 yields 
If < CY,{tm-tk+ir'^^ {s~h)ds 

k=0 •^*'= 

m— 1 

= CAt3/2^(m-fc-l)-l/^ 



k=0 



We can bound the sum above by 2M^/^, therefore we have 
/f + /f^ < C{At + Ati/2+«) < C(Ai2«). 

Let us estimate If. Using the fact that Sh is bounded and Assumption ^ yields 
(with R = PhF"[X{tk) + r[0{s) - 0{tk))){0{s) ~ 0{tk), 0(5) - O(tfe))) 



t23 



2n1/2 



< 

m-1 rtk 






fc=o ■"^fe 

< E / is~tkf ds 

< CiAtf". 

Combining if + if and if yields the foUowing estimate 

E {WhW^y^'^ < C{l\f% < C(Ar). 
Combining the previous estimates for the term / yields for Xq E D{—A), 

m-l ptk + i 



Ell/ 



fc=0 •'^'= 



|2)l/2 < dh'+Af'+Y^I (E||X(ifc)-X 

and for Xa e V{{-A)'i) 

2n1/2 



fcll'^ 



1/2, 



Ell/ll 



< C(C/2(;,2 ^ ^^.) ^ ^ / (E||X(tfe) - X 

/c=0 •^*'' 



.^IP)'^') 



,1/2 



, E||//|| 



,1/2 



Finally we combine all our estimates on /, // and /// to get (E||/| 

1 /? 
and (E||///|P) and use the discrete Gronwall inequality to complete the proof. 



3.3. Proof of Theorem 2.4 for (SETDl) scheme. We now prove convergence 

/ \ 1/2 

in II^{U) and estimate (E||X(tm) — -^mll//i(n)) ■ ^^^ the proof we follow the 

same steps as in previous section for Theorem 2.2 and we now estimate (3.8) in the 
H^ norm. 
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2.2 



in Section 



3.2 



using 



Let estimate (E||//||^wf^j)^/^. As in the proof of Theorem 

the re gularity of the noise 0{t) e V{-A) = ij2(0) nH, Vt ep, T] and the^opcrty 

(2.10) of the projection Ph yields 

E||//||2,,(,,) = E\\PhPN{Oit„,))-PN{0{t„MHHn) 

< Ch^-E\\{PN{0{UnmHHn) 

< Ch^-E\\Oitr,MH-in) 

< Ch^E\\i-A) Sit,n~s)dW{s)f 

Jo 

< Ch^ / \\i-A)S{t„,-s)\\lods 

Jo ^ 

< Ch" f \\i-A)Q'/YHsds 

Jo 

< Ch'^, 

thus (E|l//|l^,(^))i/2 < Ch. 

Using the regularity of the noise again and the equivalency ||.||hi(si) = ll(~^)^^^- 
we also have 



= E||(I-P^)(-A)-i(-A)iO(i„)||?,,(^) 

= m-Ay/^{i-p^){-A)-\-Ayo{tm)rH.^n) 

< \\i-A)'^\l-PN)i-A)-'fE\\{-A)0{Un)f 

< \\{-Ay/\l-PN){-A)-'fE\\{-A)0{Un)f 



< C[ inf A, 

We now estimate the term / from ( 3.8 ) in the H^ (il) norm noting that from ( |3.9[) we 
have / = Ii+ 12+ 13 + /4 . Estimates on /i follow immediately from equations ( 3.5 ) 
and K^ of Lemma O and then for /i, if X^ E V((-A)'-''+^^/'^) C V, equation 



(3.3) of Lemma 3.1 



gives 



and if Xo e V{{-A)) == H n H'^ifl), 

mhrmin))'^' <ch [nxowi 



1/2 



2(0) 



1/2 



If F satisfies Assumption [3] (b) , then using the Lipschitz condition, the triangle 
inequality, the fact that Ph is an bounded operator and Sh satisfies the smoothing 
property analogous to S{t) independently of h [T3], i.e. 



Shit)v\\l,,^)<Ct-'^^v\\ veVh t>0, 
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we have 

k=a •'*■'' 
'"^1 ftk+l 
< CJ2 {t„r-s)-'/HnF{X{tk))-F{{ZJi + PhPNO{tk))\\'y/^)ds 

m-l ptk + 1 

Once again using Lipschitz condition, triangle inequality, smoothing property of 
Sh, but with Lemma |3.2| gives 

(E||/3ll^.(n))'/' 

< Y. {nShit„^-s)Ph{F{x{s))-F{x{t,mj,,^^^y^')ds 

< C^E / (t™-s)-i/2(E||^(X(s))-f(X(ffe))||)i/2d5 

< Cj] / it^-s)-'/'iE\\X{s)^X{t,)\\l,^^^y/^ds 

k=0 •'*■'' 

. f E||Xo||^+i + (^E^sup^||F(X(s))||hi(o)) + 1 j 

N-l/2 






fc=0 " ^'= 
2 



. E||Xo||^+i + E sup \\F{X{s))\\H^n) + 1 

Y V 0<s<T 

As in the previous theorem, we use the fact that 

™-l rtk+i 



I — n "^ ik 



k=0 •"■>' 

Then if Xo e P((-^)('^+i)/2) we have finally found 

{nhfrn^n))"^ <C{Mr/^ [nM?,+i+ (^ sup |lF(X(.))|Ui(o)) +1 

V V 0<s<T 



1/2 



In the same way, if Xq G X>(-A) we obviously have {'E\\I:i\\]J^^^^■^f/'^ < C(At)i/2-«; 
by taking 7 = 1 — e in Lemma |3.2[ e > small enough. 
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If F{X) e V then by (3.5) of Lemma 3.1 we find 



"-l rtk+i 

/"i-l rtk + i \ / \ 1/2 

< Ch[J2 {t^-s)-'/^ds]( snj>nF{X{smmm 



.fc=0-'= / ^°^^^^ 



< Ch. 



Combining our estimates, for F{X) e F we have that : ii Xq e 2?((-A)(''+i)/2) 
then 

™-l /-ifc + i 

If Xo e X'(-A) then 

where C > depending of the T, the initial sohition Xg, the mild solution X, the 
nonlinear function F . 

Combining our estimates ( E||/||^wf^^ j , f E||//||^i.j^, j and f E||///||^l.f^^ j 
and using the discrete Gronwall lemma concludes the proof. 

3.4. Proofs for the (SETDO) scheme. Recall that 



- 5,,(i„)F,,Xo +Y.[ ^''(*™ - tk)PhF{Yt)ds 

k=o ^•^*'= 

= 5;,(i„)F;,Xo +Y.[ ^"(*™ - tk)PhF{Yt)ds ) + F;,PArO(i„) 
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As in the Theorem |2.2| we obviously have 

= Zt„ + o(i„) - y,;^ 

- Z(t„0 + Oitra) - {Zt + PhPNO{tm)) 

- (X(f„0 - Z,;;J + {PN{0{t„,)) ~ PhPN{0{t^))) + (O(t„0 - PAr(0(i„))) 

= I + 11 + III. 

(3.10) 

The proofs are therefore as [TJ Theorem 2.6 and Theorem 2.7] but with 5*"^^ 
replaced by Sh{tm — tk) and using the similar estimates as in the proofs of Theorem 



2.2 and Theorem 2.4 for the (SETDl) scheme. 



4. Implementation & numerical results 

4.1. Efficient computation of the action of (pi, i — 0,1. The key element in the 
stochastic exponential schemes is computing the matrix exponential functions, the 
so called </?;— functions. It is well known that a standard Pade approximation for a 
matrix exponential is not an efficient method for large scale problems [3U |3H1 ES] ■ 
Here we focus on the real fast Leja points and the Krylov subspace techniques 
to evaluate the action of the exponential matrix function ipi{At Ah) on a vector 
V, instead of computing the full exponential function ipi{At A^) as in a standard 
Pade approximation. The details of the real fast Leja points technique and [351 
I5il [TTj for the Krylov subspace technique are given in [SH |3I1 132] . We give a brief 
summary below. In [TT we have compared the efficiency of the two techniques for 
deterministic advection- diffusion-reaction. 

4.1.1. Krylov space subspace technique. The main idea of the Krylov subspace tech- 
nique is to approximate the action of the exponential matrix function ipi [At Ah) on a 
vector V by projection onto a small Krylov subspace Km = span |v, AhV, . . . , A^~ v| 
The approximation is formed using an orthonormal basis of V™ — [vj^jVj, . . . ,y^] 
of the Krylov subspace K^, and of its completion V,„_|_i — [V„j, v„j^]^] . The basis 
is found by Arnoldi iteration which uses stabilised Gram-Schmidt to produce a se- 
quence of vectors that span the Krylov subspace. 
Let ef be the i*^ standard basis vector of W . We approximate (pi{At Ah)v by 

(4.1) f^{AtAh)v « ||y||2V„+i^,(AtH™+i)e™+^ 

with 

H,„+i = „ „"' "I where H„ = V^,A,,V„ = [ft-ij]. 

Y U, • • • ,U, llm+l,m U J 

The coefRcient hm+i.m is recovered in the last iteration of Arnoldi's iteration [M1IT71 
I33| . For a small Krylov subspace (i.e, m is small) a standard Pade approximation 
can be used to compute ipi{At'H.m+i) , but a efficient way used in [S^ is to recover 
ipi{AtIim+i)&T~^^ directly from the Pade approximation of the exponential of a 
matrix related to Hm |34j . In our implementation we use the functions expv . m and 
phiv.m of the package Expokit |34j, which used the efficient technique specified 
above. 
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4.1.2. Real fast Leja points technique. For a given vector v, the real fast Leja points 
approximate ipi{At A}i)v by P„j(At A/j)v, where P^ is an interpolation polynomial 
of degree m of (pi at the sequence of points {Ci}"lg called spectral real fast Leja 
points. These points {Ci}"La belong to the spectral focal interval [a, f3] of the matrix 
At Ah, i.e. the focal interval of the smaller ellipse containing all the eigenvalues of 
At Afi- This spectral interval can be estimated by the well known Gershgorin circle 
theorem |40j . In has been shown that as the degree of the polynomial increases and 
hence the number of Leja points increases, convergence is achieved |39j . i.e. 



(4.2) 



lim \\ipi{AtAh)v- P,n{AtAh)y\\2 = 0, 



where ||.||2 is the standard Euclidean norm. For a real interval [a, /3], a sequence of 
real fast Leja points {^i}i=o ^^ defined recursively as follows. Given an initial point 
^0, usually ^o = /?, the sequence of fast Leja points is generated by 

-1 

(4.3) TTl6-efe|= max^TT U-^fel J = 1,2, 3,. 



J — i J- 

k—O ' k—Q 



We use the Newton's form of the interpolating polynomial Pm given by 

(4.4) Praiz) = if, Ko] + y^ip^ Ko , 6 , ' ' ' , C,] H (^ " ^^^^ 






fe=0 



where the divided differences Lpi [•] are defined recursively by 



(4.^) 



Vi[^]] =V^{£.3) 



Vi[^],i 



J,?J + 1, ■ 



,6] 



Vi Ki + 1, Cj+2, • • ■ iik]~ ^i [?j, Cj + 1, ■ • • , Cfc-l] 



An algorithm to compute the action the action of the exponential matrix function 
Lpi{At Ah) on a vector v can be found in [17] where the standard way is used to 
computer the divided differences. Due to cancellation errors this standard way 
cannot produce accurate divided differences with magnitude smaller than machine 
precision. Here we used the efficient way to computer the divided differences [39ll42j. 
In 111] it is shown that Leja points for the interval [—2, 2] assure optimal accuracy, 
thus for the spectral focal interval [a, /3] of the matrix AtA^, it is convenient to 
interpolate, by a change of variables, the function (^^(c + 7^) of the independent 
variable ^ e [-2, 2] with c = (a + ,3)/2 and 7 = (/? - a)/4. It can be shown [12] 
that the divided differences of a function /(c + 7^) of the independent variable ^ 
at the points {^i}j^Q C [—2,2] arc the first column of the matrix function /(L,„), 
where 



■L'TT?, ^ Cljn 



■l^r. 



Li,,, — 



1 



1 



v 



1 u I 



We then conclude that the divided differences of Lpi{c + 7^) of the independent 
variable i G [-2,2] at the points {6}^o ^ ["2,2] is (^,(L„,)e;"+^ where e5"+^ is 
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the first standard basis vector of M'"+^. Taylor expansion of order p with scaling 
and squaring is used in [521132] to compute (y9i(L„)e™^^- In practice the real fast 
Leja points is computed once in the interval [—2,2] and reused at each time step 
during the computation of the divided differences. We use the efficient algorithm 
of Baglama et al. [31] to compute the real fast Leja points in [—2, 2]. 

4.1.3. Numerical construction of noise. We relate the decay of the eigenvalues qi of 



Q in (1.2) to the covariance function and discuss implementation. For concreteness 



we examine A on [0, Li] x [0, L2] with Neumann boundary conditions. For the noise 
in H^ , r = 1,2 we take the following values for {qij}^, ^q in the representation 



(1.2) 



(4.6) %,, =r/(i+jr, r>0. 

We call noise in H^ when the eigenvalues satisfy ( |4.6[ ). Consider the covariance 
operator Q with the following covariance function (kernel) with strong exponential 
decay [28l|29] 



Cr{(xi,vi);{x2,y2)) = 



ibih 



exp 



{X2 



xif , {y2-yif 



h\ 



bl 



where 61 , 62 are spatial correlation lengths in x— axis and y- axis respectively and 
F > 0. This covariance function is frequently used in geosciences to generated 
a random permeability (see |17l I30j). It is well known that the eigenf unctions 
{e^ e, }i.j>Q of the operator A = DA with Neumann boundary conditions are 
given by 



.(') 



1 

Li' 



A« = 0, 



.« 



^os(Af).), Af) = f 



with I G {1, 2}, i = 1, 2, 3, • • • and corresponding eigenvalues {Xi,j}ij>o given by 



In order to put the noise W in form of the representation (1.2 1, let us give the 
following lemma. 



Lemma 4.1. Let b and A he two real numbers. We have the following statement 
exp ^7 75" I I cos(Aa::)da; = 26 exp 



Proof. Note that 

r' + OO 

exp 

1 
2 



-f (^ ) )cos(Aa;)dx 



+ 00 



( J^- 



■ X 
4P 



— 1X3 



452 



-\-iXx 



--{Xbf 

IT 



dx 



1 - 
2' 



(XbY 



+ CXJ 



/ 



26 



n Xb 

x — i — — 

\/7r 



V^ Xb 

-— 2;+j— ^ 

2b Jtt 



dx. 
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Since for any contour (C) in complex plane we have ^„ exp(— z^)(iz = 0, taking (C) 
to be a rectangle with vertexes in complex plan —a, a,a + id, —a + id yields 

/a /'d /•a /'d 

e-'^^dx + i e-^^+'y^' dy - e-^''+"^^' dx - t e-^-^'+'^^'dj/. 
-a Jo J -a Jo 

Since 

pd pd j*d 

1/ e-(±°+'^)'d?;| = I / e(-'^'±'2ay+y=)^y| < g-a^ / ^y' ^y ^ q ^j^en a -^ +00 
Jo Jo Jo 

then when a — > +00, we have 



00 



e"" dx= e-("+''^) dx for all d G 

) J —00 

Using previous results allows us to have finally 



{Xbf ^^^ _{^ 

26 

4 V 62 



exp ( ~7 ( 77 ) ) cos{Xx)dx = e tt j e \ y dx 



= 2be TT 



—00 
2 



by using the fact that /_ e ^ dx = y/n. 



D 



Recall 24 that the covariance operator Q may be defined for / £ L'^{il) by 



g/(.T)= / Criix,y)fiy)dy. 
Jn 

Assume that the eigenfunctions of the operator Q are the same as the eigenfunctions 
of —A, for bi <C Li and using the strong exponential decay of Cr we have : 

46162/ / Cri{xi,yi);{x2,y2))cos{X'l^^X2)cos{xf\2)dy2dx2 
Jo Jo 

^ /-^^ / vr f {x2~x,f \\ (1) 

= r / exp -- p I I cos(A> 'X2)dx2 



-l^e.J-l 



{y2 - yif 



cos{xf^ y 2) dy 2 



= r / exp ( -- ( -2- I ) cos(A|^^(a; + a;i))(ix 



x/ expf--f-^jj cos( A^.^^ {x + yi)dx) 

r / exp I -- f -3- j j cos(Af (x + a;i))da; 



— 00 



X / '^^P ( ~ 4 ( ]^ ) ) *^°^('^i ^^ + yi))c'2; 



46162 cos(A« xi) cos(Af yi) T exp ("-^ ((Af ^61)2 + {Xfb2)')) . 
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It is important to notice that in the previous expressions we have used the fact that 



+00 



exp 



■K I X 



(i)„ 



cos(A!^ 'x)dx 



by Lemma |4.1| and 



+ C30 



exp 



■K I X 

4 Ifof 



2b i exp 



Ai) 



{>^fhf 



ie{i,2} 



sin(A^ 'x)dx — 



because the integrand is an odd function. Then the corresponding values of {qi.j } 



in the representation ( 1.2 1 is given by 



i+i>o 



qt; 



Fexp 



During our simulation, the process 



O, 



Jtk 



-'t)A, 



dW'^T) 



is generated in Fourier space as in [7] by applying the Ito isometry in each mode, 
which yields 

(4.7) (e., O.)) = e-^-^* (^ (l - e"^^'^*)) '^' i?,.. 



i e Xjv = {1, 2, 3, ..., N} , fc = 0, 1, 2..., M — 1 and Ri^k are independent, standard 
normally distributed random variables with means and variance 1. For efficient 
computations we use the inverse fast Fourier transform or some variant : eg for 
Neumann boundary conditions we use the inverse discrete cosine transform. 
The exponential functions in the schemes (SETDO) and (SETDl) are computed 
either using the real Leja points technique or the Krylov subspace technique. For 
noise with exponential correlations, bi> 0,i — 1,2 we have \\{—AY/^Q^/^\\hs < 00, 
r ~ 1,2. Furthermore AssumptionHis obviously satisfied with V — H — L^{fl) and 
6 = 1/2. We therefore expect the higher temporal order, i.e. close to 1 with initial 
data Xq = when F is taken to be linear. We need to consider the projection P^ 
of the noise onto the computational grid. There are two cases. When the vertices 
of our finite element mesh matches the evaluation points of the noise term 0{t) the 
projection Ph is trivial. We also used the centered finite volume [TT] discretization. 
Here Ph is trivial when the center of every control volume is an evaluation point 
0{t). Of course in general the evaluations points of the noise term 0{t) do not 
necessarily need to match the finite volume or finite element grids. In this case the 
noise needs to be regular for a good projection (see assumption |4]). 
In our simulations we examined both a finite element and a finite volume discretiza- 
tion in space and take as a domain fl = [0, 1] x [0, 1]. For time discretizations we 
compare the schemes here with an semi-implicit Euler Maruyama method (denoted 
'Implicitfem') and the semi-implicit Euler Maruyama of [T] that uses linear func- 
tionals of the noise as in (4.7 1. We denote by 'Implicitfem' the graph for standard 



semi-implicit with finite element method for space discretization with exponential 
correlation function,'SETDlfem' and 'SETDOfem' the graph for schemes (SETDl) 
and (SETDO) with finite element method for space discretization with exponen- 
tial correlation function, 'Implicitfvmr', r = 1, 2 the graph for standard implicit 
with finite volume method for space discretization with H^ noise, SETDlfemr and 
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SETDOfemr, r = 1,2 the graph for the schemes (SETDl) and (SETDO) with 
finite element method for space discretization with H^ noise, SETDlfvmr and 
SETDOfvmr, r = 1,2 the graph for the schemes (SETDl) and (SETDO) with 
finite volume method for space discretization with H^ noise, Modifiedlmplicitfvmr, 
r = 1, 2 graph for the modified implicit scheme constructed in [1] with finite volume 
method for space discretization with H^ noise. 

4.1.4. A linear reaction- diffusion equation. As a simple example consider the reac- 
tion diffusion equation in the time interval [0, T] with diffusion coefficient D > Q 

(4.8) dX^{DAX -\X)dt + dW X(0)=Xo, 

with homogeneous Neumann boundary conditions in Q. Here A is a constant related 



to the reaction and in the notation of (1.1 1 F{u) — ~Xu and obviously satisfies 
condition (a) of Assumption (Is]). For this linear equation we can construct an exact 
solution up to any spectral projection error. We compute the exponential functions 
ifi with the real fast Leja points technique. The absolute tolerance used is 10^^. 
We start by examining in Figure [T] convergence with H'" noise, r = 1,2. The 
figure compares the finite element discretization for schemes (SETDO), (SETDl), 
the standard implicit Euler-Maruyama scheme and the modified implicit scheme 
introduced in [1] which also uses a linear functional of the noise. We observe that 
schemes with finite element and finite volume space discretization have the same 
order of accuracy. In Figure [l] (a) the noise is in H^ and the diffusion coefficient 
\s D = 1. We clearly see improved accuracy of the schemes that use the linear 
functions of the noise : namely (SETDO), (SETDl) and modified implicit over 
the standard semi-implicit method. Not only is there an improved constant but the 
temporal order is higher. Numerically we find from Figure [T] an order of 0.97 for 
(SETDO), (SETDl) and for the modified semi-implicit Euler-Maruyama scheme, 
which are in excellent agreement with the theoretical value of 1 from the theory, 
the order of convergence of the standard implicit scheme is 0.30. We also see that 
the scheme (SETDO) and the modified implicit scheme have approximately the 
same order of accuracy and that (SETDl) is slightly more accurate comparing the 
schemes (SETDO) and the modified semi-implicit Euler-Maruyama. In Figure fl] 
(b) the noise is H^ and diffusion coefficient D = 1/100. The error here is dominated 
by space discretization error, as a consequence to see the convergence with need 
small Aa; and Ay. We observe again that the schemes using the linear functionals 
are more accurate. We also see from both Figure [T] (a) and (b) that (SETDl) 
is slightly more accurate than (SETDO) by some constant. The temporal order 
of convergence for schemes using linear functional of the noise is 0.97 and 0.5 for 
standard semi-implicit scheme. From Figure [I] (a) to Figure fll(b) we observe that 
as the noise is regular the gap between errors in different schemes become small. 
In Figure[2]we show results with the exponential covariance function for the noise, as 
the noise is certainly in if, r = 1 or 2 we expect a rate of convergence close to one. 
The figure compares the finite element discretization for schemes (SETDO) and 
(SETDl) against the standard implicit scheme. The temporal order of convergence 
of the schemes (SETDO) is 0.80 and (SETDl) is 1.05 and 0.80 for standard implicit 
scheme. We see the improved accuracy in the schemes (SETDO) and (SETDl) 
comparing to the standard implicit. We also see the better accuracy of the scheme 
(SETDl) compared to (SETDO). 
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(a) 



(b) 



Figure 1 . Convergence in the root mean square L^ norm at T = 1 
as a function of At with H'', r = 1,2. (a) Shows convergence for 
finite element and finite volume discretizations with r = 1, D = 1, 
A = 1, r = 1 and Ax = Ay = 1/100. In (b) we show convergence 
for finite element and finite volume discretizations with r — 2, 
D = 1/100, A = 1, r = 1, Ax = Ay = 1/400 (smah to have a good 
look of convergence) . The initial data is Xq — and the simulation 
is for (|4.8|) with 20 realizations. 




Figure 2. Convergence in the root mean square L^ norm at T = 
1 as a function of Ai with exponential covariance function with 
D = 1, X = 0.5, r = 1 and regular mesh coming from rectangular 
grid with size Ax — Ay = 1/100. The simulation is for (4.8) with 
correlation lengths bi = I 
is given by Xq — 0. 



0.2 and 10 realizations. Initial data 



4.2. Stochastic advection diffusion reaction. As a more challenging example 
we consider the stochastic advection diffusion reaction SPDE 



(4.9) 



dX = DAX - V • (qX) - 



X 



X + l 



dt + dW, 



with mixed Neumann-Dirichlet boundary conditions, and constant velocity q — 



(1,0) for homogeneous medium. In terms of equation ( 1.1 ) the nonlinear term F is 
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given by 

(4.10) F{u) = -W ■ {qu) - -—-, ueR+ 

and clearly satisfies Assumption [s] (b). For heterogeneous medium we used three 
parallel high permeability streaks. This could represent for example a highly ide- 
alized fracture pattern. We obtain the Darcy velocity field q by solving the system 

r Vq = 

(4.11) <^ fc(x)„ 

with Dirichlet boundary conditions T]-, — {0, l}x [0, 1] and Neumann boundary F^y 
(0, 1) X {0, 1} such that 

1 in {0} X [0, 1] 
' = I 
and 



^ 1 in {Li} X [0, 1] 



~kS/p{x,t) -11 = in F]y 

where p is the pressure, /i is dynamical viscosity and k the permeability of the 
porous medium. We have assumed that rock and fluids are incompressible and 
sources or sinks are absent, thus the equation 

"fc(x), 



(4.12) V-q==V 



-Vp 



= 



comes from mass conservation. 

To deal with high Peclet flows we discretize in space using flnite volumes. Simula- 
tions are in L^{n) since the discrete L^{il.) norm is easy to implement for all types 
of boundary conditions. We can write the semi-discrete flnite volume method as 

(4.13) dX'^ - {AhX'' + PhFiX'^) + b{X'')) + PhPNdW, 

where here Ah is the space discretization of DA using only homogeneous Neumann 
boundary conditions and b{X'^) comes from the approximation of diffusion flux at 
the Dirichlet boundary condition size. 

We compute the exponential functions ifi with Krylov subspace technique with 
dimension m = 6 and the absolute tolerance 10~^ and the real fast Leja points 



technique for ipo- In Figure 3(a) we shows the convergence of schemes (SETDO), 
(SETDl) and standard implicit scheme with H^ noise for homogeneous medium, 
the 'true solution' is the numerical scheme with smaller time step At = 1/15360. 
All the schemes have 1/4 for temporal order of convergence. We can also observe 
the accuracy of the scheme (SETDl) and (SETDO) comparing to and standard 



implicit scheme in Figure 3(a) In Figure 4(a) we shows the convergence of schemes 



(SETDO), (SETDl) with H noise for heterogeneous medium. The two schemes 
have the same error. The corresponding mean of CPUtime for the scheme (SETDO) 
is given in Figure [4(d) | We observe a slightly efficiency gain using the Leja points 
technique compared to the Krylov subspace technique during the evaluation of the 
action of ipo ■ 

In conclusion we obtained superior convergence for the stochastic exponential in- 
tegrators using linear functionals of the noise with a flnite element discretization. 
Furthermore we have shown that these schemes that require the exponential of a 
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non-diagonal matrix can be efficiently implemented for finite element and finite 
volume discretizations of realistic porous media flow with stochastic forcing. 









^^ 






^^ Implicitfvm? 
line wrlh slope V4 


Xy^ 




/ 

^ 


^ 
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(a) 



(b) 



Figure 3. (a) Convergence of the root mean square 1? norm at 
T = 1 as a function of Ai with 30 realizations with Ax = Ay = 
1/160, Xq = 0, r = 0.01 for homogeneous medium. The noise is 
white in time and in H"^ in space, r = 2. The temporal order of 
convergence in time is 1/4 for all schemes. In (b) we plot a sample 
'true solution' for r = 2 with At — 1/15360. 
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Figure 4. (a) Convergence of the root mean square L^ norm at 
T = 1 as a function of At with 30 reahzations and Ax = Ay = 
1/160, Xq = 0, r = 0.01 for heterogeneous medium. The noise is 
white in time and in H^ in space, r = 2. The temporal order of 
convergence in time is 0.26 (close to 1/4) for the two methods. In 
(b) we plot a sample 'true solution' for r = 1 with At = 1/15360. 
In (c) we plot the streamline of the velocity field. In (d) the mean 
CPUtime for SETDO using the Krylov and Leja points. The Peclet 
number for the flow is 16.58. 
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